Flame propagation in random media 
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We introduce a phase-field model to describe the dynam- 
ics of a self-sustaining propagating combustion front within 
a medium of randomly distributed reactants. Numerical sim- 
ulations of this model show that a flame front exists for re- 
actant concentration c > c* > 0, while its vanishing at c* 
is consistent with mean-field percolation theory. For c > c* , 
we find that the interface associated with the diffuse com- 
bustion zone exhibits kinetic roughening characteristic of the 
Kardar-Parisi-Zhang equation. 
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The behavior of a nonequilibrium system is often lim- 
ited by a reaction front between one metastable or un- 
stable phase and another more stable phase. A common 
but spectacular example of this occurs in slow combus- 
tion where a flame front forms and can propagate ■ 
Here we study a phase field model of combustion in the 
absence of convection and apply methods developed in 
the study of phase transitions to the asymptotic behav- 
ior of a self-sustaining combustion front growing within 
a medium of randomly distributed reactants. Both the 
formation of the front as well as its universal dependence 
on length and time are examined. We find numerically 
that the combustion front exists for reactant concentra- 
tion c>c* >0 in <i = 2 dimensions, where the behavior 
at c* is consistent with that of mean-field percolation. 
For c > c* , we find that the diffuse combustion zone 
exhibits kinetic roughening characteristic of the Kardar- 
Parisi-Zhang interface equation 0. 

While we expect our phase-field model to describe the 
universal features of combustion or reaction fronts in the 
absence of convection we shall motivate it below in terms 
of a specific example, forest fires. The physics associ- 
ated with forest fires has recently received attention (^-^] 
due to the potential relationship with the concept of self- 
organized criticality, introduced by Bak |5] and collabo- 
rators. In most cases studied to date, cellu ar automaton 
models on a lattice have been used ■ In these works 
a collection of forests which can burn and subsequently 
reappear is considered. In contrast, this paper focuses 
on systems in which the reacting element cannot sponta- 
neously reappear. 

Our model consists of two coupled reaction-diffusion 



equations, one governing the concentration of reactants, 
and the other the dynamics of a thermal field. Unlike dis- 
crete lattice models this model incorporates the interplay 
between long-range thermal diffusion and local random 
concentration fields. Within our model, variations in the 
local temperature field T(x, t) at position x and time t, 
are due to three effects: (i) thermal diffusion through 
the medium; (ii) Newtonian cooling due to coupling to 
a heat bath; and (iii) generation of heat, limited by ac- 
tivation, from the reactants. Explicitly, the temperature 
field obeys 



dT 
~dt 



DV 2 T-T[T-T ]-V-WT + R(T,c), (1) 



where D is the diffusion coefficient, T is the thermal dissi- 
pation constant, and To is the constant background tem- 
perature of the bath responsible for Newtonian cooling. 
For completeness we have included convection due to an 
external forcing V, but we shall hereafter set this term 
to zero. 

Nonlinearities enter through the reaction rate R(T, c), 
which is limited by the local concentration of reactants 
c. Activation implies that this term is proportional 
to exp(— A/T), where A is an activation energy, and 
Boltzmann's constant has been set to unity. Further- 
more the rate is limited by the local flux cx VT, so 
that the local energy produced per unit time is q(T) — 
T 3 / 2 exp(— A/T), where the additional factor of T sets 
the scale of energy. In the combustion literature the 
multiplicative prefactor is T Q , where a — 0(1) varies 
according to the conditions present. However, the par- 
ticular form is relatively unimportant since the dominat- 
ing temperature dependence is the exponential activa- 
tion. Hence, on measuring temperature in units of the 
activation energy, we model the reaction rate as 



(2) 



where Ai is a dimensionless constant. This completes 
the formulation. Initial conditions determine the ran- 
dom distribution and concentration of reactant. Here we 
have initially distributed the reactant at random, with 
the probability that a given system site is occupied be- 
ing c. 
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For the remainder of this paper, we consider d = 2, 
where a front initially parallel to the y axis propagates 
in the x direction. The dimensionless parameters are set 
to D = 0.2, T = 0.05, T = 0.01 and Ai = 8, and time is 
measured in units of those for the reaction, A2/A1, and 
length in units of the dimension of the reactant. In our 
numerical work, the mesh size in space is set to Aa: = 1, 
while the mesh size in time is Ai = 0.01; tests of smaller 
mesh sizes give qualitatively similar results. It is useful to 
relate these choices of parameters to the specific example 
of a forest fire. For example, the constant A2 in Eq. 2, can 
be found in terms of the density and specific heat of air 
and the activation temperature of wood p0| ] . In physical 
units, we have D ~ lm 2 s -1 , T ~ 0.05s" 1 , T - 10K and 
c p ~ 5Jg~ 1 K~ 1 and A ~ 500K. With the exception of 
To, these are comparable to real systems. Our small To 
has been chosen to give enhanced cooling and hence keep 
diffusion fields relatively short ranged. This allows us to 
perform our numerical integrations with good accuracy 
without having to simulate extremely large systems. Test 
runs show that our results are relatively insensitive to the 
choice of T Q . 

Due to the activated nature of the combustion process, 
a self-sustaining propagating combustion front requires a 
sufficient amount of heat to be released during combus- 
tion. The source for this heat is the reactant concentra- 
tion c. Since activation limits the production of this heat, 
we expect existence of a nonzero concentration c* , below 
which the fire will spontaneously burn out due to insuffi- 
cient heat production. That is, for c < c* the velocity of 
the front v(c) = 0, while v is nonzero for higher concen- 
trations. For quantitative analysis the position h(x,t) of 
the front, where v = dh/dt, is defined as the position x 
where the temperature is maximum at a given time and y 
position. The variable h is then a single-valued function 
of y. 

For large concentrations, v is constant after an initial 
transient, and decreases with c. The transient increases 
as c* is approached. In the vicinity of c*, the asymptotic 
velocity approaches the relationship v(c) ~ (c — c*)*, 
where (f> is an exponent similar to that obtained in per- 
colation theory M . 

To determine the scaling exponent in the case of a 
random background, Eqs. (2) and (3) were numerically 
solved on a lattice using periodic boundary conditions in 
the y direction and fixed boundary conditions in the x 
direction. The dimension of the system is L in the y di- 
rection, and well exceeds vt max in the x direction, where 
tmax is the maximum time studied. At every site, the den- 
sity variable c is initially either zero or one, and randomly 
distributed with average c. The fire is started at the far 
left by igniting a complete row of "trees" at y = 0. Af- 
ter a short transient, the propagating fire front assumes 
a steady-state average velocity v(c). In Figs. 1(a) and 
(b) typical configurations of the propagating tempera- 
ture field are shown at c = 0.65 and 0.225. For lower 
densities, the front becomes very irregular and finally 
stops propagating. Calculating the velocity numerically 



for L = 200 we find c* = 0.19 ±0.02 and <f> = 0.46 ±0.09. 

Mean-field theory is useful to understand these perco- 
lation results. Consider a uniform distribution of reac- 
tants, whose density variable c is every where equal to a 
constant. In this description there are no longer varia- 
tions in T in the y direction. Assume there exists a mean- 
field temperature front T m moving with constant velocity 
v m (c). Using dT/dt = v m dT/dx, we obtain a nonlinear 
consistency relation between T m and v m (c). We have 
solved this mean field model numerically, revealing a de- 
pendence of v m (c) in c of the form v m (c) oc (c — c*)*, 
near c* = 0.19, where <fi — 0.5. The exponent 0=1/2 
is also a consequence of the following argument. Let %t 
be a point sufficiently ahead of the peak of T m , beyond 
which dissipation exceeds activation in the thermal dif- 
fusion equation. Expanding R(T) to linear order around 
Ttaii = T m (xt), we find that the leading edge of T m goes 
as T m ~ exp[(-v m /2D - 0^ - 4D(X lC q'(T tail ) - T)x}. 
The requirement that T m does not develop any oscilla- 
tory components gives v m > 4D\iq'(T tai i)(c — c*) 1 / 2 , 
with c* = r / \q' (T ta ii) Ai) . This analysis is analogous to 
the method used in Ref. |ll[] to find front velocities in 
the context of epidemic models. 

To incorporate finite-size effects, we use a scaling form 

w(c,X) ~X-^fi[( c -c*)Z 1/l/ ]. (3) 

This is the same as that used in percolation theory ||. 
Here v is the correlation- length exponent £ ~ (c — c*) _ly , 
and the scaling function f2(x — ► 00) ~ ar. In Fig. 2, we 
show numerical results for v(c, L)L^/ V .vs. (c — c*)L 1 / v 
for nine different system sizes ranging from L — 4 to 
L = 200. Using c* = 0.19 and <j> = 0.46, we find that the 
best collapse occurs for v = 0.6 ±0.1, as shown in Fig. 2. 

The results for the critical exponents are consistent 
with the mean field exponents of percolation, for which 
<j) = v = 1/2 ||. Qualitatively, heat propagation in our 
model is limited by a percolation lattice, provided by the 
random density field c. Below c* , the connected cluster 
available for front propagation breaks down, and the fire 
spontaneously dies out. The mean-field nature of the 
critical exponents is due to the long-range nature of the 
diffusion field associated with T. 

For c > c* , it is clear from Fig. 1 that the propagating 
interface associated with T develops large fluctuations 
and appears rough. We define the width of the interface 
by to = ((h — (/i)) 2 ) 1 / 2 . Rough interfaces often satisfy 
the scaling relation @Q w(t,L) - ^/(t/L 2 ) for large 
L and t. An important example of this is the Kardar- 
Parisi-Zhang (KPZ) interface equation 0, for which the 
exact and nontrivial exponents are (3 — 1/3, z — 3/2, and 
X — z(3 — 1/2, in d = 2. In our case, for any given value 
of c > c* , we expect the width to obey this scaling form, 
i.e., for large t, we expect w ~ in the limit t <C L z , 
where finite-size effects can be neglected. To account 
for the dependence of the width on concentration in this 
limit, we propose 

w{c,t)=^{c)w s {t/r{c)), (4) 
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where w s (t — > oo) ~ t 13 , and for c near c*, £(c) ~ 
(c — c*)~ v and r(c) ~ (c — c*)" A , where A is a slow- 
ing down exponent. In mean field, 2v = A = 1. It is 
straightforward to generalize this form to include finite- 
size effects. In Fig. 3 we show the scaled width w s plotted 
.vs. the scaled time t s — t/r, for seven different values 
of c with L = 200. For this L, finite-size effects play no 
discernible role. The inset shows the original data set. 
A transient time to nas been subtracted, which has been 
determined from the point where v(c) reaches a constant 
value. From the fitted £(c) and r(c) for the data collapse 
we cannot accurately estimate v and A, although they 
are consistent with the mean-field values. 

From the scaled data of Fig. 3, we determine the rough- 
ening exponent [3. The running slope of the data from a 
logu; vs. logf plot gives an effective /3(£), which is shown 
in Fig. 4. After an initial transient the slope clearly tends 
towards (3 — 1/3, which is the exact KPZ value. We 
have also analyzed the data by calculating the difference 
w(bt) — w(t) = A(b@ — l)tP, where & is a constant (e.g., 
b — 2). From this we find j3 = 0.34 ± 0.04, which is our 
best estimate of this exponent. 

In the limit where the width saturates due to finite-size 
effects, i.e., t » L z , the scaling relation gives w(c,L) ~ 
L x . Using system sizes L = 50, 76, 100, 150, 200, 300, 
400 and 600, we obtain \ — 0-5 ± 0.1 for c = 0.5 and 
X = 0.5 ±0.2 for c = 0.85, as compared to the exact KPZ 
value of x = 1/2- Our results for (3 and x are therefore 
in good agreement with those of the KPZ equation ^,|). 

To summarize, in this work we have developed a re- 
alistic phase-field model for combustion fronts. We find 
a percolation transition at a critical density c* w 0.19, 
below which the fire will spontaneously die. We have 
analyzed the nature of this transition and found critical 
exponents which are consistent with those of mean-field 
percolation. Furthermore, above c*, we found that the 
diffuse combustion front displays kinetic roughening. By 
an appropriate generalization of the usual scaling form 
for interfaces, we have shown that the roughening expo- 
nents are compatible with the KPZ universality class. 

The Centre for Scientific Computing (CSC) Co. of Es- 
poo, Finland has provided most of the computing re- 
sources for this work. This work has also been supported 
by the Academy of Finland, the Natural Sciences and En- 
gineering research Council of Canada and les Fonds pour 
la Formation de Chercheurs et I'Aide a la Recherche de 
Quebec. 



[3] P. Bak, K. Chen, and C. Tang, Phys. Lett. A 147, 297 
(1990). 

[4] B. Drossel and F. Schwabl, Phys. Rev. Lett. 69, 1629 
(1992); P. Grassberger and H. Kantz, J. Stat. Phys. 63, 
685 (1991). 

[5] P. Bak, C. Tang, and K. Wiesenfeld, Phys. Rev. Lett. 
59, 381 (1987). 

[6] Recently, J. Zhang, Y. C. Zhang, P. Alstrom, and M. T. 
Levinsen, [Physica A 189, 383 (1992)] have interpreted 
the results of burning sheets of paper in terms of kinetic 
roughening. They found \ ~ 0.71, which is larger than 
we find. This may be due to the correlated distribution 
of fibers within paper. 

[7] M. Kardar, G. Parisi, and Y. C. Zhang, Phys. Rev. Lett. 
56, 889 (1986); for a review see J. Krug and H. Spohn, in 
Solids Far from Equilibrium: Growth, Morphology, and 
Defects, ed. G. Godreche (Cambridge University Press, 
Cambridge 1991). 

[8] D. Stauffer, Introduction to Percolation Theory, 2nd ed. 
(London, Taylor & Francis Ltd, 1985). 

[9] D. Stauffer, Phys. Rep. 54, 3 (1979). 
[10] For a rough estimate of the parameters entering our 
model for the specific case of a forest fire, we assume trees 
burn via a steady-state reaction with air, treating air as 
an ideal gas. The flux of molecules striking the wood sur- 
face is N s ~ (n/4) y/{2T/m), where n is the number den- 
sity, and m is the mass per molecule in air. Since the num- 
ber of molecules reacting is N r = N s e~ A ^ T , the energy Q 
released per unit tree area and time equals the energy of 
the molecules entering the reaction, Q = (3T/2)N s e A ^ T . 
If r is a tree radius, the heat given off per unit volume 
of a tree site is 2Q/r. Writing T in units of A, the tem- 
perature change per unit time corresponding to this heat 
production is (2Q/r) /(c p p), where c v is the specific heat 
and p is the mass density of air. Now, on recalling that 
we are scaling temperature in units of A, we obtain the 
parameter, A 2 = (3 V / 2/4rc p )(l/m) 3/2 (A) 1/2 . 
[11] J. D. Murray, Mathematical Biology (Springer Verlag, 
1989). 

[12] F. Family and T. Vicsek, J. Phys. A 18, L57 (1985). 

FIG. 1. The temperature field T(x,y,i) for a moving fire 
front in a uniform, random forest with (a) c = 0.65, and (b) 
c = 0.225. The black pixels correspond to the temperature 
field; the higher the temperature the blacker the pixels. The 
interface h(x, i) is defined by the curve outlined by the black- 
est pixels. The light grey pixels to the right of the interface 
represent reactant. 



FIG. 2. Finite size scaling of v(c,L). The 

main figure shows ln(v(c, L)U / '' / ) .vs. ln((c — c*)L 1 ^ I/ ). 
The inset shows the unsealed data for system sizes 
L = 4,6,8,24,44,54,64,104,200, from right to left. Sizes 
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FIG. 3. Crossover scaling function w a plotted .vs. 
t s — t/r. The inset shows the concentration dependent width 
w(c,t), for c = 0.4,0.5,0.6,0.7,0.75,0.80 and 0.85 from top 
to bottom. The roughness increases with decreasing density. 
A transient time t and the corresponding offset w has been 
subtracted from each w(c, t) 

FIG. 4. A log-log plot of the scaling function w s of Fig. 
3. The inset shows the effective /3 as a function of time. The 
straight line represents f3 = 1/3. 
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